Fluid mixing optimization with reinforcement learning

Fluid mixing is crucial in various industrial processes. In this study, focusing on the characteristics that reinforcement learning (RL) is suitable for global-in-time optimization, we propose utilizing RL for fluid mixing optimization of passive scalar fields. For the two-dimensional fluid mixing problem described by the advection–diffusion equations, a trained mixer realizes an exponentially fast mixing without any prior knowledge. The stretching and folding by the trained mixer around stagnation points are essential in the optimal mixing process. Furthermore, this study introduces a physically reasonable transfer learning method of the trained mixer: reusing a mixer trained at a certain Péclet number to the mixing problem at another Péclet number. Based on the optimization results of the laminar mixing, we discuss applications of the proposed method to industrial mixing problems, including turbulent mixing.

Fluid mixing plays a fundamental role in various industrial processes. However, most mixing processes are empirically designed by using trial-and-error methods through physical experiments, rather than mathematical optimization. Although turbulence is an "effective mixer" 1 , in some cases (e.g., a bioreactor or a mixer in food industry processes), turbulent mixing is not appropriate because strong shear flows damage the materials to be mixed. Moreover, sustaining turbulent flows in micro-mixers is difficult due to low Reynolds numbers; this necessitates enhanced mixing by laminar flows. Therefore, mixing optimization by laminar flows is crucial. Several analytical studies have evaluated the efficiency of laminar mixing protocols 2-5 , e.g., proving the exponential bounds of the mixing speed; however, the research on constructive optimization methods remains limited.
This study proposes a mixing optimization based on reinforcement learning (RL) as a constructive method. To illustrate the effectiveness of the RL algorithm for fluid mixing optimization, we first summarize its mathematical framework. The RL algorithm is formulated in terms of the Markov decision process (MDP) 6,7 : M = {S, A, p 0 , P, R} , where S denotes the set of states, S = {s 1 , · · · s |S| } ; A denotes the set of actions, A = {a 1 , · · · a |A| } ; p 0 denotes the probability distribution of the initial state, p 0 : S → [0, 1] ; P denotes the transition probability, P : S × S × A → [0, 1] ; and R denotes the reward function, R : S × A → R . The initial state, s 0 , is determined by p 0 (·) , and in the next step, the state is determined by the transition probability, P(·|s 0 , a 0 ) , which requires the action, a 0 . The action is determined by the policy, π : S → A , as a = π(s) . The RL algorithm is implemented to determine the optimal policy, π * , for the given MDP, which maximizes the expectation of the cumulative reward, ∞ t=0 γ t R t+1 . Here, γ ∈ (0, 1) denotes the discount factor and R t+1 := R(s t , a t ). The RL algorithm maximizes the cumulative reward (i.e., global-in-time) rather than the instantaneous reward, R t (i.e., local-in-time). Therefore, it is suitable for global-in-time optimization problems. Designing efficient mixing protocols is one of the global-in-time optimization problems, as the final scalar field depends on the temporal order of actions in the entire mixing process, which includes stretching and folding by fluid flows and its coupling with molecular diffusion. An illustrative example was presented in History matters of Villermaux 8 . Despite the effectiveness of the RL algorithms in solving a diverse range of problems in fluid mechanics [9][10][11] , including nuclear fusion 12 and turbulence modelling 13 , the fluid mixing problem remains unexplored.
The RL algorithm is suitable for the global-in-time optimization problems, but not for problems with a high-dimensional state space in general, which is known as the curse of dimensionality 6 . Particularly, the high dimensionality of the state space for fluid mixing renders the implementation of the RL algorithm challenging. This study investigates an optimization problem formulated by Mathew et al. 2 , in which the velocity field is given by the superposition of the prescribed fields. This reduces the dimension of the state space for the fluid motion to one 2 ; a single parameter, denoted by θ later, determines the state of the fluid motion. This optimization problem was based on a physical experiment using the electromagnetically driven flow 14 . The conjugate gradient descent method was introduced as a prototype of the fluid mixing optimization 2 . To ensure that the RL algorithm can handle the flow field with a reduced degree of freedom, we focus on the same optimization problem. www.nature.com/scientificreports/ In this paper, we demonstrate for the first time that the RL algorithm is suitable for fluid mixing optimizations. This algorithm identifies an effective flow control, which results in exponentially fast mixing without prior knowledge. The mechanisms behind efficient mixing are uncovered by focusing on the flow around the fixed points from the viewpoint of dynamical system theory 15,16 . This study also proposes an effective transfer learning method for the trained mixer by considering the diffusion effect on mixing. Based on the optimization results of the laminar mixing, we discuss applications of the proposed method to industrial mixing problems, including turbulent mixing, in the "Conclusion and discussion" section.

Methods
Fluid mixing optimization. We consider the following optimization problem formulated by Mathew et al. 2 as the benchmark problem, in which the velocity field, u(x, y, t) = α 1 (t)u 1 (x, y) + α 2 (t)u 2 (x, y) , is used. Here, u 1 (x, y) = (− sin(2πx) cos(2πy), cos(2πx) sin(2πy)) and u 2 (x, y) = u 1 (x − 0.25, y − 0.25) (see Fig. 1a). The time evolution of the passive scalar, c(x, y, t), is described by the advection-diffusion equations on the twodimensional torus, T 2 (the periodic boundary condition): where Pe ∈ (0, ∞] represents the Péclet number. As a constraint on the flow control, the time integral of the kinetic energy, 1 . We set α(t) = 2 E (cos θ(t), sin θ(t)) , by which the constraint is always satisfied. We also set E = 1.25 as in Mathew et al. 2 . In this problem, the velocity field, u(x, y, t), is determined by a single parameter, θ(t) , referred to as the flow parameter.
The variance of the scalar field is often used to measure the degree of mixedness. However, as it is a conserved quantity in the absence of diffusion (i.e., d dt T 2 c p dx ≡ 0 (∀p ∈ N) ), it is unsuitable as a measure of the mixing process. Instead, we employ the mix-variance defined by �(c) = �c� 2 H −1/2 := k 1 √ 1+(2π�k�) 2 |c k | 2 , where c k denotes the Fourier coefficient of the scalar field 17 . The mix-variance is equivalent to the Mix-Norm that was originally introduced to characterize the multiscale property of the mixed scalar field 17 . Furthermore, Mathew et al. 17 showed the equivalence between the decay of �(c) , the weak convergence in L 2 , and the mixing of ergodic dynamical systems in Theorem 3.2 (see also Lin et al. 3 for the extension of the theorem). To summarize the optimization problem, we use the RL algorithm to determine the function, θ : [0, 1] → R , that minimizes the mix-variance at the end of the mixing process, �(c(·, t = 1)).
We conduct a numerical simulation of the advection-diffusion equations (Eq. 1) using the fourth-order Runge-Kutta scheme for the temporal integration with t = 0.001 and the Fourier spectral method for spatial discretization with a grid of 250 × 250 , which is the same as the one used in Mathew et al. 2 .
Optimization of mixing with RL. Here, we consider the optimization of the action-value function (Q instead of the policy π , and obtain the optimal Q func- www.nature.com/scientificreports/ tion, Q * : S × A → R . The Banach fixed-point theorem mathematically ensures that such an optimal Q function exists as a fixed point of the Bellman operator 6,7 . We obtain the optimal policy as π * (s) := argmax a∈A Q * (s, a). As a standard implementation of the RL algorithm, we employ the deep Q network 18 , which approximates the Q function by using the deep neural network denoted by Q w : R N s × A → R . Here, N s and w denote the dimension of the state space and the connection weights in the neural network, respectively. The inputs to the network are the scalar field, c(x, y, t), and the velocity field, u(x, y, t). The values of these fields on T 2 are observed on the square grid 83 × 83 , and the state, s, of the MDP is defined as the observed values of the velocity field, , and those of the scalar fields over the last five steps; that is, and t O denotes the time interval of the successive observations. Therefore, the dimension of the state space is N s = 7 × N O . The network consists of four hidden layers, and each activation function is ReLU as Mnih et al. 18 . The discount factor is γ = 0.99 . The more details of the deep Q network structure and its implementations are described in the "Supplementary Information". The initial distribution, p 0 , is given by the delta function such that θ(0) = 0 and c(x, y, 0) = sin(2πy).
The time interval of the successive observations is t O = 0.004 , which is the same value used in the benchmark problem 2 , and t Q = 5 t O , where t Q denotes the time interval of the successive updates of the Q function. Hence, for each period of t Q , the RL algorithm observes the scalar fields determined by the advection-diffusion equations (Eq. 1) with the fixed velocity field. Then, the Q function, i.e., the weights in the neural network, is updated. A single unit of episode corresponds to a single mixing process, i.e., solving the initial value problem of the advection-diffusion equations (Eq. 1) for 0 ≤ t ≤ 1 . The total number, N e , of episodes for the training is N e = 4000 . The results with the larger number of episodes, N e = 5000 , are qualitatively identical to those with N e = 4000.
As the action, A, of the MDP, the RL algorithm can change the value of the flow parameter, θ(t) (0 ≤ t ≤ 1) . The velocity field, u(x, y, t), is determined by the single-parameter θ(t) , and the flow control is realized by changing θ(t) . The discretization of the temporal change of the flow parameter is where ω + = π/(4�t Q ) and ω − = −π/(4�t Q ) . The action, ω , is selected following the ε -greedy method 6,7,18 , which changes the value of ε linearly from 1 to 0.001.
The reward function, R, is defined by using the mix-variance, , which is set to be a monotonically decreasing function of to ensure that the smaller value of represents a better mixed scalar field: where ˜ , 0 , and T denote a threshold, an initial value, and a target value of the mix-variance, respectively. By definition, R = −1 initially, and R = +1 if the mix-variance, , attains the target value. The values of ˜ and T are set based on the Péclet number:

Results
RL-based optimization. The optimization results are presented in the absence of diffusion ( Pe = ∞ ). The optimal policy, π * : R N s → A , approximated by the deep Q network, is obtained from the RL-based optimization. Thereafter, the state vector, s t ∈ R N s , determines the optimal action through ω t = π * (s t ) . This determines the velocity field during the next interval, t Q , which advects the scalar field, and the process carries on to the next observation. This flow controller based on the optimal policy, π * , is referred to as a trained mixer. Figure 1b shows from the left to right panels that the trained mixer makes the scalar field, c(x, t), evolve in time. Here, the black and white colors correspond to the high and low values of the scalar field, respectively. The trained mixer produces a complex layered structure of the scalar field. The following subsection presents a detailed description of the successive stretching and folding of the interface between the two colors.
The mix-variance, � n (t) (n = 1, . . . , 4000) , is shown in Fig. 1c. During the initial stage of the training, (i.e., in the first half of the total episodes such as n = 1, 800, and 1600), the RL algorithm with the ǫ-greedy method chooses actions randomly. Although this "random mixer" can decrease the mix-variance, such mixing is inefficient, as explained below.
Mathew et al. 2 reported that the proposed conjugate gradient descent method resulted in �(t = 1) ≃ 6 × 10 −3 ; this value of the mix-variance is used for the comparison as a reference. In the first half of the total episodes, the mix-variance at the end of the mixing process, � n (t = 1) , is larger than the reference value; that is, the insufficient training of the mixer results in inefficient mixing. Conversely, � n (t = 1) is reduced in the latter half of the total episodes, n = 2400 , 3200, and 4000. Particularly, 3 × 10 −3 < � n (t = 1) < 4 × 10 −3 for n = 4000 , which are almost identical to (slightly smaller than) the reference value. Interestingly, the mix-variance decreases exponentially fast for 0.3 ≤ t ≤ 1 for the latter episodes such as n = 3200 and n = 4000 . While here we focus on the quantitative comparison using the mix-variance, there are some qualitative differences between the method by Mathew et al. 2 and our RL-based method. In the "Conclusion and discussion" section, we illustrate significant advantages of the RL-based method. Figure 1d presents the mix-variance at the end of each mixing process, � n (t = 1) , which fluctuates due to the ε-greedy methods and the fact that the policy, Q w , is not converged. However, the fluctuation decreases as the episode progresses; see also Figs. S1 and S2 in the "Supplementary Information". The RL algorithm significantly decreases the mix-variance, � n (t = 1) ; that is, the RL-based optimization effectively enhances the mixing.
Characteristics of the trained mixer. The flow parameter in episode n is denoted by θ n (t) . In the first half of the training, n < 2000 , the flow parameter, θ n (t) , evolves randomly in time due to the ε-greedy methods www.nature.com/scientificreports/ and the fact that the policy is not converged. However, as the episode progresses, θ n (t) converges to a single function, θ * (t) , except for the final stage of the process, as shown in Fig. 2a. The time series of θ n (t) consists of square waves, as the velocity field (i.e., θ n (t) ) is fixed in each interval, t Q . The optimal mixing process by the trained mixer corresponding to θ * (t) is divided into the following three stages:  www.nature.com/scientificreports/ flow is temporally periodic in the middle stages, which are shown in the lower panels in Fig. 2b. Each velocity field has eight fixed (stagnation) points, u 1 and u 2 . Half of them are elliptic; that is, the Jacobian matrix has purely imaginary eigenvalues. The other half are saddle points; that is, the Jacobian matrix has real eigenvalues 15,16 . We focus one of them at (x, y) = (0.5, 0.5) , which is depicted by the red point in each panel of Fig. 2b as a reference. The material line around the fixed point is stretched along the unstable eigendirections when the fixed point is a saddle, whereas it is folded (approximately π/2 rotation) when the fixed point is elliptic. The local stretching and folding around the eight fixed points occur simultaneously, resulting in efficient mixing. The use of the specific protocol by the trained mixer with the constant angular frequency, θ(t) = ω * t , is explained in the section of "Conclusion and discussion" section.
Remarkably, the period of the flow in the middle stage, 2π/ω * , which determines the period of the successive switching of the saddle and elliptic types of the fixed points, is optimal in the following sense. Apart from the RL algorithm, we conduct numerical simulations of the scalar field advected by the flow determined by θ(t) = ωt with a constant angular frequency, ω , throughout the mixing process, 0 ≤ t ≤ 1 . The inset of Fig. 2a shows �(t = 1) evaluated for ω ∈ [0, 30] . The minimum of �(t = 1) in this setting is obtained at ω ≃ ω * . This implies that the RL algorithm determines the optimal angular frequency, ω * , without any prior knowledge, and the trained mixer uses the temporally periodic flow with the optimal period in the middle stage of the process.
Comparison with randomized mixers. To characterize the flow by the trained mixer in the initial and middle stages, we introduce three different mixing processes, called randomized mixers: • Completely randomized mixer: It uses the random controller that takes one of the three actions, ω ∈ {0, ω + , ω − } , independently, with the same probabilities for all the stages ( 0 ≤ t ≤ 1). • Partially randomized mixer I: It uses the trained mixer for the initial stage ( 0 ≤ t < 0.3 ), and then switches to use the random controller for 0.3 ≤ t ≤ 1. • Partially randomized mixer II: It uses the trained mixer for the initial and middle stages ( 0 ≤ t < 0.7 ), and then switches to use the random controller for 0.7 ≤ t ≤ 1.
Numerical simulations are conducted 200 times independently for each control. Fig. 2c presents the probability density functions (PDFs) of the mix-variance, �(t = 1) , at the end of the mixing process. The gray solid line indicates the value of the mix-variance by the trained mixer, � n (t = 1) (n = 4000) (see Fig. S1 in the "Supplementary Information" for the related PDF of the trained mixer).
The top panel of Fig. 2c depicts the PDF in the case of the completely randomized mixer, where the mixvariances are larger than the reference value of the trained mixer. The left and right panels of Fig. 2d represent the final state of the scalar field, c(x, t = 1) , produced by the trained mixer and one completely randomized mixer that presents the mix-variance, �(t = 1) , close to the median value of the PDF. Videos 1 and 2 in the "Supplementary Information" correspond to the scalar fields mixed by the trained mixer and the completely randomized mixer, respectively. Large unmixed blobs remain in the scalar field produced by the completely randomized mixer. That is, the training mixer with the RL algorithm is effective. The second panel of Fig. 2c depicts the PDF in the case of the partially randomized mixer I, which is more effective than the completely randomized mixer. However, a substantial gap exists between the results of the partially randomized mixer I and that of the trained mixer. This indicates that the mixing process during the middle stage is also crucial. Finally, the third panel of Fig. 2c depicts the PDF produced by the partially randomized mixer II. The results are almost identical to those obtained using the trained mixer. Therefore, the effectiveness of the partially randomized mixer II is the same as that of the trained mixer. These observations demonstrate that the mixing process during the initial and middle stages is essential for the mixing efficiency, whereas the mixing process during the final stage is not.
Transferring the trained mixer. This subsection considers the diffusion effect on the RL-optimization of the mixing described by the advection-diffusion equations (Eq. 1) with finite Péclet numbers. The details of the problem settings are identical to those in the previous sections, except for the values of the Péclet numbers. The RL-based optimization is applied to the mixing problem for the case of Pe = 10 2 , 10 3 , and 10 4 , which are as effective as for the case of Pe = ∞ , regardless of the Péclet numbers. For instance, at Pe = 100 , the mix-variance, � n (t) , decreases faster for the later episodes, as shown in the inset of Fig. 3b, where n = 1, 600, 1200, 1800, 2400 , and 3000 and lighter (thicker) curves correspond to larger n. We note that the curves of � n (t) for n ≥ 1200 are almost the same, implying that the RL algorithm converges to find the optimal policy at n = 1200 . Interestingly, this convergence is faster than the case of Pe = ∞ (Fig. 1c). The number of episodes required for convergence is n ≃ 3000 at Pe = ∞ ; however, n ≃ 1200 seems to be sufficient for convergence around Pe = 100.
The diffusion effect appears in the flow controls at the later stages. If the mixer successfully generates fine layered structures in an early stage, the flow control becomes less important in the later stages of mixing due to the diffusion effect. In other words, at a low Péclet number, once the RL algorithm finds the optimal mixing control in an early stage of mixing, nothing is to be learned as the diffusion reduces the mix-variance rapidly, irrespective of the control by the mixer. This may result in the faster convergence observed above. The implications of the fast convergence at the low Péclet numbers to the training mixer are given in the "Conclusion and discussion" section.
This diffusion effect implies the asymmetric transferability of a trained mixer; that is, a mixer trained at a high Péclet number can be used for mixing at a lower Péclet number, whereas the converse does not hold true. Let Pe T be the Péclet number where the mixer is trained, and the asymmetric transferability is then rephrased as follows: the trained mixer can be reused for the same mixing process for the range of (0, Pe T ] . Figure 3a  www.nature.com/scientificreports/ and the thin red lines indicate the results for the case of Pe T = 100 . In Fig. 3a,b, the solid, dashed, and dasheddotted lines indicate the results with different random numbers for learning. The mixers that trained at Pe T = ∞ realize the exponentially fast mixing for the entire process when we use it for Pe = ∞ . On the other hand, the mixers that trained at Pe T = 100 realize the exponentially fast mixing only for the first half of the process, but fail to mix for the latter half. Fig. 3b presents the mix-variance, �(t) , for 0 ≤ t ≤ 1 at Pe = 100 . Similar to Fig. 3a, the thick blue lines represent the results for the case of Pe T = ∞ , and the thin red lines represent the results for the case of Pe T = 100 . Unlike the case of Pe = ∞ , no significant difference exists between the results for the cases of Pe T = 100 and Pe T = ∞ , and both cases realize the exponentially fast mixing. In summary, the mixers of Pe T = ∞ can be used for the mixing at Pe = 100 , whereas the converse does not hold true. Therefore, a mixer trained at a higher Péclet number can be used for the mixing process for a broader range of Pe.

Conclusion and discussion
By illustrating why the RL algorithm is suitable for fluid mixing optimization, we demonstrated, as a proofof-concept, that the mixer trained by using the RL algorithm is effective for the two-dimensional fluid mixing problem (Fig. 1), which paves the way for the development of RL-based training of mixers. The proposed method was quantitatively evaluated by focusing on the benchmark problem of the mixing optimization studied in the pioneer work 2 . In addition to the comparison of mix-variance values, we note that our RL-based method solves the optimization problem in more restrictive conditions compared to the method proposed by Mathew et al. 2 . For example, in our setting, the number of the states of the velocity field is restricted to eight, θ = 0, π/4, . . . , 7π/4 . Furthermore, the proposed method is more flexible; that is, it uses only the scalar and velocity field as the input to the neural network. Provided that these fields can be observed, physical implementations are possible in principle, even if the evolution equations of these fields are unknown. For instance, mixing problems of granular or viscoelastic fluids are essential; however, the evolution equation of such a complex material has not necessarily been established, and therefore, the conjugate gradient descent method 2 cannot be applied to these industrially fundamental problems. On the other hand, the RL based method is equation-free, hence it is applicable if the sensory data of the mixture states are available as the input to the neural network.
The optimized mixing process was divided into three distinct stages. It is particularly interesting to note that, in the middle stage, the optimized flow is temporally periodic with the constant angular frequency. Here, we discuss why the RL algorithm makes the angular frequency constant. The fixed points in both velocity fields, u 1 and u 2 , are located at the same position and are homogeneously placed in the domain, T 2 . If the angular frequency is not constant, the switching period between the saddle and elliptic types of the fixed point may differ at each location. This spatial difference makes the scalar field inhomogeneous. The inhomogeneity increases the amplitude of the Fourier coefficient of the small wavenumber, thereby increasing the mix-variance. Consequently, the time variation of the angular frequency results in the larger value of the mix-variance. The RL algorithm employs the constant angular frequency to avoid this undesired effect. The detailed justification of the aforementioned interpretation is one of the future works.
Another related future work is to understand the optimal mixing in more detail. For instance, we claim that the random variation of the flow parameter in the final stage ( t > 0.7 ) is not essential for optimal mixing, in the sense that the results of the partially randomized mixer II (Fig. 2c) and the trained mixer (Fig. S1 in the "Supplementary Information") are almost identical. However, there is a small difference between these PDFs, which suggests that the randomization of the actions in the final stage may eliminate some actions, which the RL algorithm regards to be essential, in the optimized mixing process.
For practical application, reduction of learning costs is crucial. Despite the effectiveness of transfer learning in reducing learning cost, its application to the problems in fluid mechanics remains limited 19 . In this regard, www.nature.com/scientificreports/ this study has introduced the physically-reasonable notion of the asymmetric transferability of the trained mixer. The demonstration in this study (Fig. 3) indicates that, in terms of transfer learning, the Péclet number of the source domain Pe T should be as high as possible, if the trained mixer is required to reuse for the broader range. If the mixer is trained at a high Péclet number, it can learn how to mix the scalar field to create the fine striped structures. If the trained mixer is transferred to a lower Péclet number, it makes the fine structures at the beginning of the mixing process. Then, smoothing such structures by diffusion reduces the mix-variance, irrespective of the actions of the trained mixer in the later stage. Therefore, the transfer of the trained mixer from a high Péclet number to a lower one is effective. Concerning another aspect of learning costs, we have found that the learning of mixing at a lower Péclet number converges faster (Inset of Fig. 3b). Therefore, if fast learning at a Péclet number is required, the Péclet number of the source domain Pe T should be as low as possible. Considering together with the discussion in the previous paragraph, the above discussions suggest a trade-off between broad transferability and fast learning; in other words, there is an optimal Péclet number of the source domain that balances these two advantages in each application. Although this study is restricted to the transfer of the trained mixer over the different Péclet numbers, future developments of transfer learning methods of trained mixers may be significant.
Large gaps exist between the mathematical toy problem discussed in this study and the existing mixing problems in the industrial processes. However, the results of this study indicate some directions overcoming these gaps. First, we discuss the implications of this study to turbulent mixing. Turbulence comprises multi-scale counter-rotating pairs of coherent vortices 20 , and strong turbulent mixing stems from the effective mixing around such pairs of vortices at each scale 1 . As observed in the transfer learning method, scalar mixing occurs from larger to smaller scales. Since the time scale of turbulent mixing is shorter for smaller scales, the total mixing efficiency is determined by the mixing at the largest scale. Thus, measuring the velocity and scalar field at the largest scale may be sufficient for the proposed training method. Despite the significant gap between laminar and turbulent mixing, the insights from the present study will be useful for training mixers with turbulent flows.
In addition, in industry, the multiphase and/or thermal flows with chemical reactions may have to be considered, which increases the complexity of the flow dynamics. In such cases, embedding prior knowledge, such as the evolution equations or some physical constraints into the RL-based optimization may be effective, as discussed in Brunton 11 . As another future task for the RL-based optimization in industrial mixing problems, it will be important to study the robustness of the mixing control with the obtained policy with respect to changes in the initial scalar field. Furthermore, whereas the deep Q network is employed as the first step in this study, a more specific and state-of-the-art implementation of the RL algorithm would be necessary for such complex flows. Extending the proposed method to incorporate knowledge on fluid mechanics and suitable RL implementation techniques can further enhance mixing even in industrial processes with laminar and turbulent flows.

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.